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We study the effects of uniform time delays on the extreme fluctuations in stochastic synchroniza¬ 
tion and coordination problems with linear couplings in complex networks. We obtain the average 
size of the fluctuations at the nodes from the behavior of the underlying modes of the network. We 
then obtain the scaling behavior of the extreme fluctuations with system size, as well as the distri¬ 
bution of the extremes on complex networks, and compare them to those on regular one-dimensional 
lattices. For large complex networks, when the delay is not too close to the critical one, fluctuations 
at the nodes effectively decouple, and the limit distributions converge to the Fisher-Tippett-Gumbel 
density. In contrast, fluctuations in low-dimensional spatial graphs are strongly correlated, and the 
limit distribution of the extremes is the Airy density. Finally, we also explore the effects of nonlinear 
couplings on the stability and on the extremes of the synchronization landscapes. 

PACS numbers: 89.75.He, 05.40.-a, 89.20.Ff 


I. INTRODUCTION 

Synchronization and coordination involve a system of coupled, autonomously interacting units or agents attempting 
to achieve a common goal Synchronization of a system emerges from the cumulative efforts of the individual 

entities, each regulating themselves based on the information they can gather from their neighbors on the system’s 
local state. The difficulties in synchronization or coordination problems are often compounded by stochastic effects 
and time delays [MU , preventing global coordination or consensus. Time delays between the state of the system and 
the reaction to that information (due to e.g. transmission, cognition, or execution) can pose significant challenges. 
Critical aspects of the underlying theory of delays have been long established in the context of macro-economic cycles 
as far back as 1935 In such cases, the description of the complex system can be reduced to a single stochastic 

variable (l^ - [l^ . Recent interest in the application of time delays to complex networks [3, in El provides fresh 
insights extending these older results. Understanding the dynamics across a complex network offers the possibility 
to optimize synchronization d , including weighted graphs d . Synchronization and coordination with 

delays has been studied in the stock market [ll , ecological systems |27l - l30j| , population dynamics [sil - fsl , postural 
sway and balance and the human brain [9l-ETI.l3l. It is also irimortant to understand critical functions of 

autonomous artificial systeins, such as congestion control in networks d S113) I39l - l4ll| , massively parallel [H, El and 
distributed computing E3) IH , and vehicular traffic [H EHH • The aim of this paper is to explore the effects of 
noise and delays on the dynamics in complex and random networks pol - fsl . specifically on the extreme fluctuations. 
Extreme fluctuations can have critical implications in synchronization, coordination, or load balancing problems, since 
large-scale or global system failures are often triggered by extreme events occurring on an individual node |54l - l56l| . 
In order to show the implications of the general theoretical results, we will cover the implications for typologically 
distinct networks. 

The scaling behavior of extreme fluctuations in the case of zero time delay has been investigated previously for 
small-world (SW) and scale-free (SF) networks [i^, as well as low-dimensional regular topologies [i^. Despite 

having more complex interaction topologies, coordination and synchronization phenomena of the former systems (as 
far as critical behavior is concerned) actually tend to be simpler than those of their low-dimensional regular-topology 
counterparts. This is because fluctuations of the relevant field-variables at the nodes are weakly correlated in complex 
networks [13, [H, . Hence, standard extreme-value limit theorems apply to the statistics of the extremes (as well as 

to those of the system-averaged fluctuations, i.e., the width) In contrast, fluctuations in one-dimensional regular 
lattices are strongly correlated, and the applicability of traditional extreme-value limit theorems immediately break 
down [ 4 ^, (as well as limit theorems for the sum of local variables 0). 
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While extreme-value theory for the scaling properties and universal limit-distributions of uncorrelated (or weakly- 
correlated) random variables is well established lO-IO . only a few results are available on statistical properties of 
the extremes of strongly correlated variables |57l. l58l|. Majumdar and Comtet obtained the distribution of extreme 
fluctuations in a correlated stochastic one-dimensional landscape only recently (with no time delays). In coupled 
interacting systems with no delays, possible divergences of the width and the extremes are associated with the 
small-eigenvalue behavior of the Laplacian spectrum (e.g., with long-wavelength modes in low-dimensional systems 
or low connectivity in complex networks) [4 !tI l55l . |57H59j. In the presence of time delays, however, singularities and 
instabilities can also be governed by the largest eigenvalues when the system is close to the synchronizability threshold 
[EiQ- To that end, we investigate finite-size effects and the universality class of the extreme fluctuations in complex 
networks stressed by time delays. 

In parallel to the sum of a large number of uncorrelated (or weakly-correlated) short-tailed random variables 
approaching a Gaussian distribution (governed by the central-limit theorem), the largest of these (suitably-scaled) 
variables converges to the Fisher-Tippett-Gumbel (FTG) [dll - f^ (cumulative) distribution, 


^max(i) - e ^ , (1) 

where x = {xmax — o,N)/bN is the scaled extreme [d^, with mean (x) = 7 (7 = 0.577... being the Euler constant) 
and variance (t| = (i^) — (x)^ = tt^/G. The expected largest value of the original variables (e.g., for exponential-like 
tails 0 ) scales as 


(a::max) = un + hn't 



( 2 ) 


Note that corrections to this scaling are of (!I(l/lniV), which can be noticeable in finite-size networks that are 
computationally feasible in our investigations. For comparison with numerical data, it is often convenient to employ 
the extreme-value limit distribution of the variable scaled to zero mean and unit variance, y = (xmax ~ {^max))/ 


max 


( 2 /) = e ® 


(ay+T') 


( 3 ) 


where a = The corresponding FTG density then becomes 


Pmax(2/) = ae 




( 4 ) 


We hypothesize that the FTG limit distributions of the extreme fluctuations in stochastic network synchronization 
will also be applicable to the case of nonzero time delays, provided that the large but finite system is in the synchroniz- 
able regime. Although the fluctuations at the nodes will, of course, depend on the delay, the system can be considered 
as a collection of a large number of weakly-correlated components. In contrast, in the case of a one-dimensional regular 
lattice (ring) with delayed coupling, we expect that the limit distribution of the extreme fluctuations approaches the 
Airy distribution [4511571 l58L |65| . 

Thus, provided that the system is synchronizable, the scaling with the system size and the shape and class of 
the respective extreme-value limit distributions will be the same as those of a network without time delays. To 
put it simply, time delays will impact the “prefactors” (within the syncronizable regime), but not the extreme-value 
universality class. The focus of this paper is to test the above hypotheses. 


II. EIGENMODE DECOMPOSITION, FLUCTUATIONS, AND THE WIDTH 

In the simplest linear synchronization or coordination problem in networks with delay, the relevant (scalar) variable 
at each node evolves according to 

N N 

dthi{t) = - '^Cijlhiit - r) - hj{t - t)] -f r]i{t) = -'^Tijhj{t - r) -f r]i{t) , (5) 

1=1 i=i 

where Cij is the coupling matrix and Fy = —Cij -\-6ij Cu is the (symmetric) network Laplacian. Here, we consider 
unweighted graphs, thus Cij is just the adjacency matrix and C'it = ki is the degree of node i. The noise ??i(t) is 
Gaussian with zero mean and correlations {'r]i{t)r]j (t')) = 2D6ijS{t — t')). In our simulations, without loss of generality, 
we set D = 1. We have previously studied the behavior of the average size of the fluctuations about the mean for a 
network with noise and time delays [ 1 , Qi he., the width 

(u;2(t)) = /lf][h.(t)-h(t)]A , (6) 


2 



where h(^t) = N~^ hi{t) is the mean at time t and (•) indicates averaging over different realizations of the noise. 
In the present paper, we are interested in the extremes of the fluctuations in the system at a given time. Because of 
the symmetry about the mean of the relaxation term in Eq. ([5|), the distribution of extreme fluctuations above and 
below the mean are identical, so we will reduce the presentation of results to those of the maximum of a snapshot, 
given by 


^max (0 — maX'[Aj(t)} 
i 


(7) 


where Ai(t) = hi{t) — h(t) is the fluctuation about the mean of an individual node. 

Diagonalizing T from Eq. ([5]) gives N independent modes hk{t), each of which obey an equation of the form 

dthk = -\khk{t-T)+f]k{t) , ( 8 ) 

where \k is the corresponding eigenvalue for mode k. Organizing the labels of the modes such that 0 < A/c < A^+i, 
a network with positive, symmetric couplings and a single connected component has a single (and uniform) mode 
associated with Aq = 0 , which does not contribute to fluctuations about the mean and so does not impact either 
the width or the extremes, as both are measured from the mean. The condition for the average fluctuations (/i^) to 
remain finite in the steady-state for Eq. dS]) is known exactly [1, [H, [13, [13 ll^l , 

AfcT < 7r/2 . (9) 

Hence, for the network to remain synchronizable, the above must hold for all /c > 0, or equivalently, 

T Tc — ■7r/2Aniax • (1^) 


This condition guarantees that the system avoids delay-induced instabilities and that both the width and the extremes 
will have a finite steady-state value. Further, for the simple stochastic differential equation with delay in Eq. dH]), the 
steady-state variance of the corresponding stochastic variable is also known exactly 0, 


(hi) 


nrf(X,r) ^ 

Afercos(AfeT) AfeCOs(AfeT) 


( 11 ) 


Hence, given the eigenvalues of the Laplacian for a given network, one has an exact expression for the average 
steady-state width as well is, 


N N-1 N-1 


2=1 






( 12 ) 


Of course, for a typical large complex network one does not have the eigenvalues explicitly in hand. Nevertheless, 
one can obtain them through numerical diagonalization. Hence, employing Eq. m provides an alternative to direct 
simulations of the coupled stochastic differential equations with delay Eq. ©• Equation (HID, after Taylor expansion 
of T/(Afer) in the variable r in Eq. dill) , also allows one to obtain the approximate behavior of the steady-state width 
(within the synchronizable regime t < Tc), 


7~> /-) N — 1 pv N—1 —1 2 

k—1 k—1 ^ k—1 k—1 


= {w‘^)r=0 + + O(t^) {w‘^)r=0 + Dt + 


(13) 


for large networks (l/A<Cl). In obtaining the above expression we exploited that the trace is invariant under basis 
transformation, hence, J2k=i — N{k) for an unweighted graph. The first term above is 

the width for the network with no delay, which depends strongly on the detailed structure of the graph through its 
Laplacian spectrum, {w‘^)r=o = % TEk=i ^k^ Hi- first-order correction in a network with delays is completely 
independent of any structural characteristics of the network. The second-order correction only depends on the average 
degree (average connectivity), but is independent of the local connectivity of the nodes, or the specific shape and 
heterogeneity of the degree distribution. The behavior of the width as a function of the delay for a Barabasi-Albert 
(BA) scale-free network and an Erdds-Renyi (ER) graph is shown in Fig.[l] indicating an approximately linear 
behavior for a significant portion of the synchronizable regime, in accordance with the above prediction [Eq. (|13p ]. 
For comparison, the analogous behavior of the largest fluctuations is also shown, indicating (as expected) that the 
steady-state width (ic^) and the extreme fluctuations will diverge at the same critical delay [Eg. [TU]. 
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FIG. 1: (Color online) The typical behavior of the steady-state average width and the expected extreme as a function of the 
delay r for an ER and a BA network with A^=1000 and (k) « 6. 


III. EXTREME FLUCTUATIONS 


With an understanding of the typical fluctuations of the underlying modes, we may now proceed to consider the 
extreme fluctuations in a network. Consider the covariance matrix of fluctuations at the nodes (i.e., the steady-state 
equal-time correlations), ct? = (A^Aj), and that of the modes, = (hkhe) = Ske^l (where the single subscript 
denotes the diagonal elements). The distribution of a single mode follows a zero-mean normal (Gaussian) distribution 
with a variance given by Eq. (ED, 

CTfe = (hi) = DrfiXkT) . (14) 

In turn, the fluctuations from the modes translate back to those at the nodes according to cr^ = Sar'^S~^, where S is an 
orthogonal matrix with columns composed of the normalized eigenvectors of the network Laplacian (i.e., T = SAS~^, 
where A is a diagonal matrix of the eigenvalues). Since cr^ is diagonal, this transformation can be written simply as 

{Al)=al (15) 

k 


The marginal distributions of the fluctuations at the nodes (x = Ai) are Gaussian, 


P^(x) = 




(16) 


with zero mean and variance af [Eq. dTSl) ]. The effects of several delays on the spread of the distributions for a few 
representative degree classes are shown in Fig. [51 Each panel shows the distributions for a distinct delay r, which can 
be expressed in terms of the fraction relative to the critical delay q = t/tc. 

For zero or small delays the size of the fluctuations (the width of the distributions) at a node decreases monotonically 
with the node’s degree, i.e., the larger the degree the narrower the distribution [Fig. |^a,b)]. For a sufficiently large 
delay, however, the trend changes, and the node with the largest degree can exhibit the largest fluctuations [FigE 
This can be understood by noting that in a mean-field sense, the effective coupling at each node is its degree ki 
Thus, the fluctuations at each node are approximately proportional to /(fc^r) (where the scaling function is known 
exactly [Eq. dill) ]), and it is non-monotonic in its argument SB- This trend is also illustrated in Fig. |3l For zero 
(or small) delay the average fluctuations at a node decay as {Af) ^ 1/ki UM- In contrast, the average size of 
the fluctuations as a function of the degree becomes non-monotonic for large delays [Fig. [3]. The fluctuations at the 
low-degree nodes remain largely unaffected, while fluctuations at the large-degree nodes increase significantly. 
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FIG. 2: (Color online) Individual node distributions of representative degrees in a BA network with A^=1000 and (k) « 6 for 
fractions q of the critical delay of (a) q = 0, (b) q = 0.90, and (c) q = 0.99 {q=T/ tc). 


For networks with no delays it has been established that the nodes with the smallest degree typically contribute most 
to the extremes which is still valid in the case of small delays (t/tc <C 1). For scale-free (SF) networks with 

power-law degree distributions, such as BA networks [s^, the low-degree nodes can still dominate the distribution of 
extreme fluctuations at higher delays (but r < Tc) since they are more numerous, even though the typical fluctuations 
for the highest degree node are larger than for a single low degree node. So long as the highest degree node’s 
fluctuations do not dominate the extremes of the network, the large set of lowest degree nodes will lead to the familiar 
FTG distribution for the network’s extremes [Fig. lUJa)]. Note that the approach to the FTG limit distribution can 
be very slow due to the slowly vanishing corrections for Gaussian-like individual variables Further, for larger 
delays, the convergence to the FTG density may not be monotonic due to the larger effect of the largest-degree node 
for small system sizes [Fig. |31[b)]. 

Note that the largest eigenvalue of the network Laplacian varies among individual realizations of a random network 
ensemble. Therefore, to simulate “similar” synchronization dynamics in a network random ensemble (e.g., of 1000 
networks of size Af), we kept g, the fraction of the delay relative to the critical delay, fixed in the individual network 



Further, also note that the largest 
hence it diverges with the network 
and logarithmically for ER and SW 


realizations (i.e., an individual network realization £ has a delay = 
eigenvalue of the Laplacian diverges with the largest degree in a graph 
size N in coinplex networks, e.g., in a power-law fashion for SF networks 
networks [1, Q- 

We found similar behavior for other prototypical networks [SW, ER and BA], namely for t<Tc the scaled distri¬ 
butions of the extreme fluctuations converge to the ETG density [Fig. [5]. Also, our results for ER and BA networks 
indicate that the extreme fluctuations (Amax) asymptotically approach a logarithmic scaling with the system size N 
[Fig. [ 6 ], consistent with being governed by the FTG density. 

To provide further insights to the source of the extreme fluctuations in BA networks, we also analyzed the distribu¬ 
tion of the extremes within each degree class. We have already seen that the average size of fluctuations are the largest 
for small degrees, except for near-critical delays [Fig. [3]. When the delay approaches the critical value for a given 
graph, the average size of the fluctuations becomes a non-monotonic function of the degree (in a mean-field sense, this 
behavior is naturally related to U-shape scaling behavior of the fluctuations with the effective coupling i0)- In fact, 
there is a regime where the fluctuations at the largest degree node are finite and are the largest in the network. In 
parallel with the above observation, sufficiently below the critical delay of a given graph, the extreme fluctuations will 
almost always originate in the class of nodes with the smallest degree: not only the average size of the fluctuations 
is the largest here, but also this degree class has the largest population (as given by the degree distribution). In this 
regime, it is thus expected that the global distribution of the extreme fluctuations will essentially overlap with the 
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FIG. 3: (Color online) (a) Average size of the fluctuations at the nodes as a function of their degree for a BA network 
with A^^IOOO and (fc)«6. Open symbols correspond to results based on exact numerical diagonalization of the Laplacian and 
employing Eqs. (Ha and GSl). Plus symbols (of matching colors) correspond to the direct numerical integration of the stochastic 
delay-differential Eq. m with At=5 X 10 ®. The connecting lines are the average of the degree class from these numerical 
integrations, (b) Same data as in (a) but on smaller vertical scales. 


extremes within the class of the minimum degree. Our simulations confirm this in Fig. [7Da) and (b). Further, as the 
delay approaches its critical value in the given graph, the (single) node with the largest degree will often give rise to 
the largest fluctuations in the network. This is demonstrated in Fig.[7Kc), showing that the tail of the global extremes 
coincides with the (Gaussian) fluctuations at a (single) node with the largest degree. 

When the delay in a given network is not too close to the critical delay, one can assume that the fluctuations at 
the nodes decouple. This assumption works fairly well for complex networks with no delays [1^ [13,113 • (Note that 
such assumption is ill-fated for low-dimensional spatial graphs where a large correlation length governs the scaling.) 
Now we test this hypothesis for complex networks with delays, and predict the extreme limit distribution of the 
fluctuations. The cumulative distribution of the fluctuations at a particular node i (with x = Ai) can be expressed 
in terms of the error function. 





p-x'^l2a^ 

dx'- -^ 

ai'/2TT 


1 

2 




1 

2 


11 -I- erf 



(17) 


Assuming independence of the nodes in determining the extremal fluctuations, the cumulative distribution for the 
maximum fluctuation during a given snapshot is 


P< 

max 


N 

(x) ~ 

2=1 


(18) 


The corresponding density from Eqs. m and m is then 


N 


Pmax 




(x) ~ ^P*(a;) J|F’/(a:) = 

i 


21/2-N 


N 

■E 

2=1 


-x^/2a'l 


n 


1 -I- erf 


7 V 2 , 


(19) 


The results based on the above approximation (together with those given by the actual numerical simulations) are 
shown in Fig. [S] for the distribution, and in Fig. [5] for the average of the extremes. Note that the validity of this 
approximation can break down when the delay is sufficiently close to the critical delay so that fluctuation at the few 
highest degree nodes can completely dominate the extremes. 

Also, note that the fluctuations of the mode associated with the largest eigenvalue assumes large oscillatory com¬ 
ponents for AniaxT > l/c 1,0. This manifests with the greatest amplitude at the largest degree node with strong 
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y = ( A - <A ) ) / a, 

mCLX max ^max 


FIG. 4: (Color online) Finite-size behavior of the distribution of the extreme fluctuations for BA networks with (fc)«6 and 
relative delays (a) q=0.5 and (b) 5=0.9 obtained by numerically integrating Eq. ([5]) using At=5 x 10“®. The dotted curves in 
both panels correspond to the scaled FTG density, Eq. Id- 


oscillatory components. Once the delay is sufficiently close to the critical delay so that these large-amplitude oscil¬ 
lations dominate the extremes, an additional feature emerges in the distribution of Amax. Here, there is a broader 
non-universal tail, which originates from the finite but very large fluctuations at the node with the largest degree 
[Fig. [ZKc)]. The suppression of the contribution from low degree nodes is compounded if more than one node has the 
network’s maximum degree. Periodically, when the oscillatory behavior of this node brings it back near the mean 
h{t), the global extremes can still be dominated by the FTG-distributed extremes among the lowest-degree nodes. 


IV. THE WIDTH AND THE EXTREMES ON REGULAR LATTICES 

Finally, it is worthwhile to contrast the steady-state scaling behavior of the extremes and their distributions in 
complex networks to those on regular lattices. For regular d-dimensional lattices the largest eigenvalue of the Laplacian 
is independent of the system size. Thus, for a fixed delay r, if the system is synchronizable for a particular system 
size, it is synchronizable for all system sizes, AmaxT < 7r/2, and the usual V —>■ oo thermodynamic limit can be 
considered with no delay-induced instability. Further, as iV —>■ oo, the arbitrarily small eigenvalues of the Laplacian 
(Amin ^ will dominate the sum in Eq. (HU, just like they do in systems with no delay. Hence, one can expect 

that the scaling of the width, the extremes and their asymptotic limit distribution in the synchronizable regime will 
be identical to those with no delay, governed by a diverging correlations and long-wavelength modes (associated with 
the arbitrarily small eigenvalues of the Laplacian). For illustration, we considered one-dimensional lattices (with 
nearest-neighbor coupling) with delays [Fig. |5]. For completeness, we studied the detailed finite-size behavior of both 
the steady-state width distribution P{nP) and the distribution of the extremes P(Amax)- The results of the numerical 
integration of the systems with delays (but within the synchronizable regime, Amax^ < 'x ji) show that the asymptotic 
limit distributions of the width and the extreme fluctuations approach those of the delay-free systems, the FORWZ 
distribution and the Airy distribution [13, [l^, respectively [Fig. [5]. 
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FIG. 5: (Color online) (a) Extreme fluctuation distributions and (b) rescaling of the same for several delays for SW, ER, and 
BA networks with A^=1000 and (fc}«6. The various curves correspond to the same delay fraction q for both (a) and (b). Lines 
are predictions based on exact numerical diagonalization and employing Eqs. ([Till, di), and HSl). Symbols correspond to the 
numerical integration of the stochastic delay-differential equations Eq. m with At = 5 X 10 ®. The dotted curves in all panels 
correspond to the scaled ETG density, Eq. 0. 


V. EXPLORING THE EFFECTS OF NONLINEAR COUPLINGS 

While in our current work the focus has been to understand the effects of time delays on stochastic systems with 
liner couplings, we performed some explorations on the nonlinear effects. We considered the stochastic equation 

N N 

dthi{t) =- hj) + v ^ CijCik{hj - hi){hk - hi) + rii[t) , (20) 

i=l j,k=l{j<k) 

where the effects of time delays are captured, as before, by replacing by {hi{t — T)}fLi in the right-hand side 

of the above stochastic differential equation. Such nonlinear terms can be motivated by, for example, coarse-graining 
local growth processes in networks |^. I4di - l4a[7ll| (e.g., where only local network-neighborhood minima can progress). 
Note that in one dimension, the above stochastic equation reduces to 


dthi{t) = -{2hi - h^+i - hi_i) -|- v{hi+i - hi){hi-i - hi) + r]i{t) . 


(21) 
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FIG. 6: (Color online) Scaling of the expected maximum fluctuations (Amax) with system size N for ensembles of (a) ER and 
(b) BA networks (1000 realizations). Open symbols correspond to simulating the stochastic delay-differential equations Eq. ([5]) 
with At=0.001) Solid curves show predictions based on exact numerical diagonalization and employing Eqs. (1141) . (1151) . and 

(Hi). 


The above equation is just a (somewhat unconventional) discretization of the well-known KPZ equation [zllzi, 

dth{x,t) — — iy(yh)^ + r]i(t) , (22) 

which can be easily seen by taking the naive continuum limit in Eq. dan, hi±i{t) = h{x, t)±dh/dx + ..and keeping 
only the leading-order derivatives. 

First, we studied one-dimensional regular lattices (with nearest-neighbor connections and periodic boundary con¬ 
ditions). In what follows, we parameterized the delay relative to the critical delay of the linear system for reference, 
q = T Numerically integrating the time-discretized version of the stochastic differential equation Eq. m (see Sup¬ 
plemental Material for more details), we have found that for sufficiently small values of the nonlinear coupling v and 
time delays, the system reaches a steady state with a finite width for finite systems. The width distribution [Fig. HU] 
and the distribution of the extremes for both above [A^ax = hmax — h, Fig.[TT] and below the mean [Amin = h — hmin, 
Fig. [T2| approach the FORWZ and the Airy distribution, respectively, similar to the case of pure linear couplings. 
Further, the average width and the extremes (both above and below the mean height) scale as with the system 
size [Fig. [TU]. 

Investigating the behavior of the width for larger values of the nonlinear coupling, it is clear the synchronization 
profile can diverge (the system becomes unstable), even for zero time delay (see Supplemental Material). We checked 
and tested that this instability is not an artifact of the finite time difference At, but rather it is the result of the 
non-linear term on discrete lattices in Eq. (ED). Indeed, it has been well documented that even conventional 

lattice discretization schemes of the KPZ nonlinearity can give rise to instability in a noisy environment (even though 
its spatial continuum limit is stable and the nonlinear term, in fact, exactly cancels by symmetry in the stationary 
state [t^ - I^ ). Our observed behavior, induced by the nonlinear KPZ term in Eq. ED, is just another example for 
such instability, intrinsic on discrete structures, such as lattices. 

In networks, we found similar behavior. For example, in ER networks, we found that the nonlinear term in Eq. EOj) 
gives rise to a diverging width even for our smallest value of the nonlinear coupling strength v, even in the absence of 
time delays (see Supplemental Material). We again checked and confirmed that the lack of stability in the presence of 
nonlinear couplings is not the result of insufficient time discretization of the stochastic differential equation (1201) (see 
Supplemental Material). For small values of the nonlinear coupling v, there exist, however, a long “quasi-stationary” 
period before fluctuations begin to diverge, where we analyzed the statistical properties of the extremes. We have 
found that during this quasi-stationary period, the scaled distributions of the extremes are reasonably well-described 
by the FTG limit densities [Fig. [14]. 

For SW networks, while the fluctuations (and the width) are smaller than those than on one-dimensional lattices (at 
least during the quasi-stationary period), the synchronization landscape eventually becomes unstable for sufficiently 
strong nonlinear coupling, with or without delays (see Supplemental Material). We also observe that the average 
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FIG. 7: (Color online) Contribution of a few selected degree classes for a BA network (including fcmin=3 and A:max=116) with 
N = 1000 and {k) « 6 for (a) q = 0.0, (b) q — 0.9 and (c) q = 0.99. Results are obtained from numerical integration of the 
stochastic delay-differential equations Eq. ((SJ) with At = 5 x 10“®. The numbers in brackets in the legends indicate the number 
of nodes in the corresponding degree class. 


(a) (b) 



FIG. 8: (Color online) (a) Extreme fluctuation distributions and (b) scaled version of the same data for several delays for 
a regular one-dimensional lattice with A=1000. Symbols are the results found by numerically integrating Eq. m using 
At = 5 X 10“®. The solid line corresponds to the predicted asymptotic Airy limit distribution [s^, 


width in the quasi-stationary period is decreasing with increasing average degree (fc), but at the same time, the 
duration of the quasi-stationary period is decreasing with increasing (k) (see Supplementary Material). Nevertheless, 
in the stationary (or quasi-stationary) state, the scaled distributions of the extremes are well-described by the FTG 
limit densities [Fig. [Tb]. 

It is clear that we have only begun to scratch the surface of the complexity of the behavior as a result of nonlinear 
couplings in the stochastic delay differential equations in networks. Among the important questions that one shall 
investigate are the effects of the strength of nonlinear coupling, time delay, and network size on the length of the 
quasi-stationary period. It is also clear from our explorations that uncontrolled expansions of local growth processes 
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FIG. 9: (Color online) (a) Finite-size behavior of the scaled width distribution for one-dimensional lattices for fixed time delay 
T = 0.25. The solid line corresponds to the predicted asymptotic FORWZ limit distribution [^. (b) Finite-size behavior of 
the scaled extreme fluctuation distribution for regular one-dimensional lattices for the same system sizes and delay as in (a). 
The solid line corresponds to the predicted asymptotic Airy limit distribution [^. [^. (c) Finite-size behavior of the average 
width for the same system sizes and delay as in (a). The dashed line is to guide the eyes, corresponding to the asymptotic 
theoretical prediction, (w^) ~ N. (d) Finite-size behavior of the extreme fluctuations for the same system sizes and delay as in 
(b). The dashed line is to guide the eyes, corresponding to the asymptotic theoretical prediction, (Amax) ~ 


(i.e., naive coarse graining) may result in nonlinear terms in the resulting stochastic differential equations (with or 
without delays) which give rise to instability and diverging width in the synchronization landscape. This instability 
is “real” (i.e., not an artifact of time discretization) as far as the numerical integration of the stochastic differential 
equation is concerned, but may not be present in the actual physical systems with the original microscopic rules 


l21|,l43H4a. 


VI. SUMMARY 

We have demonstrated that the extreme fluctuations in stochastic coordination or synchronization problems with 
time delays (with short-tailed node-level noise and within the linearized approximation) can fall in two main classes. 
In complex or random networks (e.g., ER, SF, or SW graphs), if the system is sufficiently large but synchronizable 
(AmaxT < 7r/2), the distribution of the extremes is governed by the FTG distribution, while the average size of 
the largest fluctuations does not grow faster than logarithmic. This type of scaling behavior can be understood as 
the fluctuations at the nodes are only weakly correlated, hence traditional extreme-value limit theorems apply. In 
contrast, in spatial graphs, fluctuations at the nodes are strongly correlated. As demonstrated on one-dimensional 
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FIG. 10: (Color online) Finite-size behavior of the scaled width distribution for one-dimensional lattices with nonlinear coupling 
V [Eq. (1211) ] and relative time delay q = t/tc- (a) u — 0.0, q = 0.0; (b) v — 0.1, q = 0.0; (c) v — 0.0, q = 0.1; (d) v — 0.1, 
q = 0.1. The solid line corresponds to the predicted asymptotic FORWZ limit distribution [^ . 


regular rings, the distribution of the extremes approaches the Airy limit distribution, while the average size of the 
largest fluctuations will scale as the width itself, e.g., as a power law in one dimension. 

Finally, we have performed some explorations on the effects on nonlinear couplings in the stochastic delay differential 
equations in networks. Our results indicate the generalized KPZ nonlinearity in discrete structures (networks) can 
ultimately give rise to instability. Even in that case, however, during a quasi-stationary period (before the width 
diverges), the statistics of the extreme fluctuations are well-described by the Airy and FTG densities, in one dimension 
and in random ER/SW graphs, respectively. Clearly, future investigations are needed to precisely characterize the 
stability conditions in the presence of nonlinear couplings in networks with (and without) time delays. 
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FIG. 11: (Color online) Finite-size behavior of the scaled distribution of the extreme fluctuation above the mean for regular 
one-dimensional lattices with nonlinear coupling v [Eq. (1211) 1 and relative time delay q = t/tc. (a) u — 0.0, q — 0.0; (b) v = 0.1, 

B : 0.0; (c) V = 0.0, q — 0.1; (d) i/ = 0.1, q = 0.1. The solid line corresponds to the predicted asymptotic Airy limit distribution 
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FIG. 12: (Color online) Finite-size behavior of the scaled distribution of the extreme fluctuation below the mean for regular 
one-dimensional lattices with nonlinear coupling v [Eq. HSU] and relative time delay q = r (a) v = 0.0, q = 0.0; (b) v = 0.1, 

g : 0.0; (c) V = 0.0, q = 0.1; (d) v = 0.1, q = 0.1. The solid line corresponds to the predicted asymptotic Airy limit distribution 

,[5i. 
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FIG. 13: (Color online) Scaling of the average width w = ^ (w^) and the expected extreme fluctuations Amax and Amin (the 
expected largest fluctuations above and below the mean, respectively) with the system size in regular one-dimensional lattices 
with nonlinear coupling v [Eq. (I21II ] and relative time delay q = t/tc- The dashed lines are to guide the eyes, corresponding to 
the scaling ~ . Note the log-log scales, (a) u — 0.0, q = 0.0; (b) u = 0.1, q = 0.0; (c) v = 0.0, q — 0.1; (d) v = 0.1, q = 0.1. 
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FIG. 14: (Color online) Extreme fluctuation distribution (in the quasi-stationary period) in ER networks with N = 100, (k) ~ 6, 
nonlinear coupling i/ = 0.1 [Eq. (1201) 1. and various relative time delays q — t (a) Distributions of the extreme fluctuations 
above the mean, Amax = hmax — h\ (b) Scaled distributions of the extreme fluctuations above the mean (with zero mean and 
unit variance); (c) Distributions of the extreme fluctuations below the mean, Amin = h — hmin; (d) Scaled distributions of the 
extreme fluctuations below the mean (with zero mean and unit variance). The solid line in (b) and (d) corresponds to the FTG 
limit distribution [Eq. (P|)]. 
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FIG. 15: (Color online) Extreme fluctuation distribution (in the quasi-stationary period) in SW networks with N — 100, 
(k) ~ 6, nonlinear coupling u = Q.l [Eq. (I20II ]. and various relative time delays q = r( tc- (a) Distributions of the extreme 
fluctuations above the mean, Amax = hmax — h; (b) Scaled distributions of the extreme fluctuations above the mean (with 
zero mean and unit variance); (c) Distributions of the extreme fluctuations below the mean, Amin = h — hmin; (d) Scaled 
distributions of the extreme fluctuations below the mean (with zero mean and unit variance). The solid line in (b) and (d) 
corresponds to the FTG limit distribution [Eq. Q]. 
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Supplemental Material 


In this Supplemental Material, we provide some additional results (figures) obtained by the numerical integration 
of the stochastic delay differential equation 

N N 

dthi{t) =- hj) + ly ^ CijCik{hj - hi){hk - hi) + r]i{t) , (SI) 

j=l j,k=l{j<k) 

discussed in Sec. V. in the main text. The effects of time delays are captured by replacing {/ii(t)}(li by — 
in the right-hand side of Eq. (jSlI) . In the above continuous-time stochastic differential equation, the noise r]i(t) 
is Gaussian with zero mean and correlations {'rii{t)rij(t')) = 2Ddij6(t — t')). In the present work, we considered 
unweighted graphs, i.e., Cij is just the adjacency matrix. The time-discretized version of Eq. (ISII) becomes [1-5] 

N N 

/ii(< -I- At) -/ii(t) =(/li -/ij) -I- Atz/ ^ CijCik{hj - hi){hk - hi) + fji{t)V2DAt , (S2) 

j=l j,k=l{i<k) 

where fji{t) are independent and identically distributed random variables for all i and t with Gaussian distribution 
of zero mean and unit variance. In our simulations, without loss of generality, we set D=l. For example, in one 
dimension, with nearest neighbor couplings and periodic boundary conditions, the above equation becomes 

hi{t + At) — hi{t) = —At{2hi — hi+i — /ii_i) -I- At v(hi+i — hi){hi-i — hi) + fii(t)V2DAt . (S3) 

Employing Eqs. (IS2I1 and (IS3I1 . we explore and display the width of the synchronization landscape as a function of time 
for a one-dimensional lattice (with nearest-neighbor coupling and periodic boundary conditions) [Fig. IS1| . ER graphs 
[Fig.[S4], and SW networks [Figs. lSTlISS] . We also demonstrate that the instability observed in the resulting stochastic 
landscapes is not the result of insufficient time discretization in the numerical integration scheme, but rather, it is 
intrinsic to the generalized KPZ coupling on discrete structures (such as lattices or networks) for sufficiently strong 
nonlinear coupling and/or large time delays. To that end, we show the average time to desynchronization (the onset of 
instability) as a function of At used in the above numerical integration scheme for a one-dimensional lattice [Fig. IS2| . 
ER networks [Fig. [S^, and SW networks [Fig. [S9]. Finally, we show the average time to desynchronization as a 
function of the nonlinear coupling ly and relative time delay q = rjTc in a one-dimensional lattice [Fig. [S3], ER 
networks [Fig. [S6], and SW networks [Fig. ISIO] . Our results for SW networks [Figs. [S7| and [S8] also suggest that 
while increasing connectivity (average degree (/c)) reduces the width in the quasi-stationary state, the quasi-stationary 
period is becoming shorter, i.e., the system is more prone to instability. 
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FIG. SI: Graphical table of the evolution of the width over time, in a one-dimensional regular lattice (with periodic boundary 
conditions) with A^=100 for KPZ-like coupling [Eq. mi Each tile shows the width (tn^) as a function of time, using identical 
scales in each tile. Data was obtained by numerically integrating Eq. (fS3l) with At = 10 Golors represent 10 distinct 
realizations of noise. Each row corresponds to the indicated value of time delay q = t/tc, and each column corresponds to the 
indicated nonlinear coupling strength v. A high-resolution version of this figure is provided separately in TIFF format among 
the Supplemental Materials (W2vsTimeMap_Ring.tif). 
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FIG. S2: Average time to desynchronization for KPZ-like coupling [Eq. (IS3ll ] at various At time steps of integration on a one¬ 
dimensional regular lattice (with periodic boundary conditions) with A=100, v = 0.8, q = 0.8. Error bars represent standard 
deviation, sampled over 100 realizations. 



V 

FIG. S3: Time to reach desynchronization for KPZ-like coupling [Eq. IMl)], as a function of the delay q = t/tc and nonlinear 
coupling strength zz, on a one-dimensional regular lattice (with periodic boundary conditions) with A=100. Data was obtained 
by numerically integrating Eq. (IS3I) with At = 10 ®. 
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FIG. S4: Graphical table of the evolution of the width over time, in ER networks with A'^^lOO and (k) « 6 using generalized 
KPZ coupling [Eq. (IS2|l ]. Each tile shows the width (w^) as a function of time, using identical scales in each tile. Data was 
obtained by numerically integrating Eq. (IS2II with At = 10“'*. Colors represent 10 distinct realizations of noise. Each row 
corresponds to the indicated value of time delay q = t/tc, and each column corresponds to the indicated nonlinear coupling 
strength zz. A high-resolution version of this hgure is provided separately in TIFE format among the Supplemental Materials 
(W2vsTimeMapER_k6.tif). 


23 


































































































































































































































































































































































































































































2 . 0 - 


1.5- 



0.5- 


0.0 . . i— . . ^ — . . ^ — . .■ 

10® lO'* 10® 10’^ 

At 

FIG. S5: Average time to desynchronization for generalized KPZ coupling [Eq. (IS2I) 1 at various At time steps of integration in 
a ER network with A=100, (fc) ~ 6, = 0.2, q = 0.2. Error bars represent standard deviation, sampled over 100 realizations. 



V 

FIG. S6: Time to reach desynchronization using generalized KPZ coupling [Eq. (EH)], as a function of the delay q = r/r^ and 
nonlinear coupling strength v, for an ER network with A=100 and (fc) « 6. Data was obtained by numerically integrating 
Eq. (IS2I) with At = 10“®. 
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FIG. S7: Graphical table of the evolution of the width over time, in SW networks with A^=100 and (k) « 3, using generalized 
KPZ coupling [Eq. (IS2|l ]. Each tile shows the width (w^) as a function of time, using identical scales in each tile. Data was 
obtained by numerically integrating Eq. (IS2II with At = 10“^. Colors represent 10 distinct realizations of noise. Each row 
corresponds to the indicated value of time delay q = t/tc, and each column corresponds to the indicated nonlinear coupling 
strength zz. A high-resolution version of this hgure is provided separately in TIFF format among the Supplemental Materials 
(W2vsTimeMapSW_k3.tif). 
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FIG. S8: Graphical table of the evolution of the width over time, in SW networks with A^=100 and (k) « 6, using generalized 
KPZ coupling [Eq. (IS2|l ]. Each tile shows the width (w^) as a function of time, using identical scales in each tile. Data was 
obtained by numerically integrating Eq. (IS2II with At = 10“^. Colors represent 10 distinct realizations of noise. Each row 
corresponds to the indicated value of time delay q = t/tc, and each column corresponds to the indicated nonlinear coupling 
strength zz. A high-resolution version of this hgure is provided separately in TIFF format among the Supplemental Materials 
(W2vsTimeMapSW_k6.tif). 
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FIG. S9: Average time to desynchronization for generalized KPZ coupling [Eq. (IS2I) 1 at various At time steps of integration in 
a SW network with A=100, (k) k, v = 0.2, q = 0.2. Error bars represent standard deviation, sampled over 100 realizations. 



V 

FIG. SIO: Time to reach desynchronization using generalized KPZ coupling [Eq. (IS2I) ]. as a function of the delay q = t jrc. 
and nonlinear coupling strength v, for a SW network with A=100 and (k) ~ 6. Data was obtained by numerically integrating 
Eq. (IS2I) with At = 10“®. 
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